In-class Exercise 5

A short description of the post.

Installing and Loading the R package

• maptools: to create spatial objects
• sf: to import geosptial data
• raster: to create raster object
• spatstat: to perform spatial point pattern analysis (SPPA)
• tmap: to plot thematic maps
• tidyverse for data handling

packages = c('maptools', 'sf', 'raster', 'spatstat', 'tmap', 'tidyverse')

for (p in packages){
if(!require(p, character.only = T)){
  install.packages(p)
}
  library(p,character.only = T)
}

Importing Geospatial Data

• Import shapefile using st_read() of sf package. Output object is in tibble sf object class.

mpsz_sf <- st_read(dsn = 'data/shapefile', layer = 'MP14_SUBZONE_WEB_PL')

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\wlwong2018\IS415 Blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

• Projection is in SVY21.

Import aspatial data from rds folder

• read_rds() of readr package is used instead of readRDS() of base R is used. This is because output of read_rds() is in tibble object.

childcare <- read_rds("data/rds/childcare.rds")
CHAS <- read_rds("data/rds/CHAS.rds")

• Note that there are some data issue in childcare data frame because Lat and Lng should be in numeric data type. The coordinate fields seem to be in decimal degree. Hence wgs referencing system is assumed.

Converting the aspatial data frame into sf objects

CHAS_sf <- st_as_sf(CHAS,
                    coords = c("X_COORDINATE", "Y_COORDINATE"),
                    crs = 3414)

• Note: st_as_sf() accept coordinates in character data type.

• lat long in decimal degrees

childcare_sf <- st_as_sf(childcare,
                    coords = c("Lng", "Lat"),
                    crs = 4326)  %>%
  st_transform(crs = 3414)

Plot to review

• alpha: to set transparency

tmap_mode("view")

tm_shape(childcare_sf) +
  tm_dots(alpha = 0.4,
          col = "blue",
          size = 0.05) +
tm_shape(CHAS_sf) +
  tm_dots(alpha = 0.4,
          col = "red",
          size = 0.05)

Geospatial Data Wrangling

Convert sf to Spatial* classes

• as_Spatial() of sf package

childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)

Convert Spatial* dataframe into Spatial* objects

• as.SpatialPoint() of as.SpatialPolygon() of maptools package

childcare_sp <- as(childcare, "SpatialPoints") # data will be dropped. 
CHAS_sp <- as(CHAS, "SpatialPoints") # data will be dropped. 
mpsz_sp <- as(mpsz, "SpatialPolygons") # data will be dropped. 

Convert from Spatial* objects into ppp object

• to drop project info. only X, Y coordinates will be left.
• as.ppp() of maptools package

childcare_ppp <- as(childcare_sp, "ppp")
CHAS_ppp <- as(CHAS_sp, "ppp")

Remove duplicated points using jitter

• to avoid points overlapping

childcare_ppp_jit <- rjitter(childcare_ppp,
                             retry = TRUE,
                             nsim = 1,
                             drop = TRUE)

any(duplicated(childcare_ppp_jit))

[1] FALSE

CHAS_ppp_jit <- rjitter(CHAS_ppp,
                             retry = TRUE,
                             nsim = 1,
                             drop = TRUE)

any(duplicated(CHAS_ppp_jit))

[1] FALSE

Note:

• must be sf or sp to plot a tmap.
• For now, we can only plot using sf plot()

Extract Punggol Planning Area

• Remember to add the comma at the end

pg <- mpsz[mpsz@data$PLN_AREA_N == "PUNGGOL",]

Convert SpatialPolygonsDataFrame into SpatialPolygons object

pg_sp <- as(pg, "SpatialPolygons")

COnvert SpatialPolygons into owin object

• owin: an object needed by spatstat to define irregular window

pg_owin <- as(pg_sp, "owin")

Extract spatial points within owin

childcare_pg <- childcare_ppp_jit[pg_owin]
CHAS_pg <- CHAS_ppp_jit[pg_owin]

• Plotting with childcare_ppp_jit would not work. Does not show the boundary of Punggol

plot(childcare_pg)

L-function

# L_childcare <- envelope(childcare_pg,
#                         Lest,
#                         nsim=99,
#                         rank=1,
#                         global = TRUE)

Code chunk for plotting interactive L-function